      PROGRAM AUTOFLAT3
!
! Experimentalni verze 3.
! predchozi verze asi moc nefungovaly!!!
!
! Toto je nejdrsnejsi verze. Pocita primo median a tak
! potrebuje opravdu hodne pameti -- typicky desitky Mb.
!

! Vyvoreni autoflat-fieldu z mnoha CCD snimku.
!
! FH, Rijen 1997
! FH, duben 1998
! FH, cerven 1998
! FH, zari 1998
! FH, rijen 1998, 7ho - vypada to velmi dobre, i drsne testovaci snimky 
!                   to uspesne projelo
!
! Porad je to jen experimenatlni verze!
! Opatrnost nutna!
!
    
       IMPLICIT NONE

       character*80 version
       parameter( version = 'Autoflat v. 0.0.3' )

       INTEGER MAXX,MAXY,MAXOBR,NSTEP
! MAXX a MAXY jsou dimenze obrazku, MAXOBR je maxiamlni pocet obrazku
! NSTEP je krok site, ve ktere se pocita median
       PARAMETER( MAXX = 511, MAXY = 341, MAXOBR = 100 , NSTEP = 1)

! SATURACE je uroven saturace
       REAL saturace
       PARAMETER( saturace = 65565.0 )

! LEVEL je vysledna predpokladana uroven medianu, kvuli zaokrouhlovacim
! chybam by nemela byt mensi nez 10 000
       REAL level
       PARAMETER( level = 10000.0 )

       CHARACTER*80 SB,BUF
       LOGICAL SIMPLE,EXTEND,ANYF,EXIST
       INTEGER I,J,L,N
       INTEGER STATUS,BITPIX,NAXIS,NAXES(2) ! dimenze snimku
       INTEGER IMAGE(MAXX,MAXY),IBUF(MAXX*MAXY)

       REAL IMAGEBUF(MAXX,MAXY,MAXOBR),fbuf(MAXX*MAXY)
       REAL IMGPAR(MAXOBR,2),medmed,medsig,r,rl,mainsig

       INTEGER BLOCKSIZE,PCOUNT,GCOUNT
       REAL BSCALE,BZERO,x,smed

       real DATUM, SKYMOD, SKYSIG, SKYSKW

       smed = 0.0

       write(*,*) version

! je aktualni pocet obrazku
       N = 0
! nacteni jmena souboru, ze STDIN
   10  READ(*,'(A)',END=20) SB   
         if( sb .eq. ' ' ) goto 10 
         if( n .eq. MAXOBR ) goto 20
         print *,sb(:50)
         n = n + 1

! otevreni souboru
         CALL FTOPEN(25,SB,0,BLOCKSIZE,STATUS)
! cteni parametru
         CALL FTGHPR(25,2,
     &        SIMPLE,BITPIX,NAXIS,NAXES,PCOUNT,GCOUNT,EXTEND,STATUS)
         CALL FTGKYE(25,'BSCALE',BSCALE,BUF,STATUS)
         CALL FTGKYE(25,'BZERO',BZERO,BUF,STATUS)
! cteni obrazku
         CALL FTG2DJ(25,1,0,MAXX,NAXES(1),NAXES(2),IMAGE,anyf,STATUS)
! zavreni obrazku
         CALL FTCLOS(25,STATUS)
! osetreni chyb pri cteni
         IF( STATUS.NE.0 )THEN
           PRINT *,'STATUS41=',STATUS
           STOP
         ENDIF    

! z kazdeho nstep-ho pixelu obrazku se udela lokalni median
        l = 0
        DO j = 1,Naxes(2),nstep
          DO i = 1,Naxes(1),nstep
            L = L + 1
            fbuf(L) = IMAGE(I,J)
          ENDDO
        ENDDO 
        call QUICK (fbuf, l, ibuf)
        CALL MMM (fbuf, l, saturace, x, DATUM, SKYMOD, SKYSIG, SKYSKW)

        IMGPAR(n,1) = skymod
        IMGPAR(n,2) = skysig

        WRITE(*,*) '       Median=',SKYMOD,l,skysig
        smed = smed + SKYMOD

! konverzni faktor intenzit na 'level' uroven
!      r = level/SKYMOD

      DO J = 1,NAXES(2)
        DO I = 1,NAXES(1)   
! prekopirovani do buferu s obrazky
          imagebuf(i,j,n) = image(i,j)
        ENDDO
      ENDDO

      GOTO 10
   20 CONTINUE 

! vypocet medianu SKYSIGu
      do i = 1,n
        fbuf(i) = imgpar(i,2)
      enddo
      l = n
      call quick (fbuf, l, ibuf)
      CALL MMM (fbuf, l, saturace, x, DATUM, medsig, SKYSIG, SKYSKW)

! vypocet medianu SKYMODu
      do i = 1,n
        fbuf(i) = imgpar(i,1)
      enddo
      l = n
      call quick (fbuf, l, ibuf)
      CALL MMM (fbuf, l, saturace, x, DATUM, medmed, SKYSIG, SKYSKW)

! prumerna hodnota SKYSIG na urovni 'level'
      mainsig = level/medmed*medsig
      write(*,*) 'mainsig=',mainsig, medmed, medsig 

! prepocet vsech obrazku tak, jako by mely stejne SKYSIG
      do l = 1,n
! intenzity se nasobi timto cislem, aby se dostala intenzita 'level'
        rl = level/imgpar(l,1)
! prepocet 'SKISIG' pro tento obrazek na standardni
        r = mainsig / imgpar(l,2)
        
        write(*,*) l,': ',r,rl 
        do j = 1, naxes(2)
          do i = 1, naxes(1)
            x = imagebuf(i,j,l) 
            x = r*(x -  imgpar(l,1)) + level
            imagebuf(i,j,l) = x
           enddo
        enddo
      enddo

c       print *,i,j,fflat(i,j)
! konecny vypocet flat-fieldu
! konstatna 'smed' udava prumernou hodnotu medianu
       smed = smed / n

       if( n .le. 1 ) stop 'Too few files to flat-fielding.'

       DO J = 1,NAXES(2)
         DO I = 1,NAXES(1)
           do l = 1,n
             fbuf(l) = imagebuf(i,j,l)
           enddo
           l = n
           call QUICK (fbuf, l, ibuf)
           CALL MMM (fbuf, l,saturace, x, DATUM, SKYMOD, SKYSIG, SKYSKW)
!           if( l .le. 10 ) print *,'l=',l
           image(i,j) = int(skymod + 0.5)
         ENDDO
       ENDDO

! cvicny vypocet medianu
        l = 0
        DO j = 1,Naxes(2),nstep
          DO i = 1,Naxes(1),nstep
            L = L + 1
            fbuf(L) = IMAGE(I,J)
          ENDDO
        ENDDO 
        call QUICK (fbuf, l, ibuf)
        CALL MMM (fbuf, l, saturace, x, DATUM, SKYMOD, SKYSIG, SKYSKW)
        WRITE(*,*) '     Konecny median=',SKYMOD,skysig,' smed=',smed,n

       INQUIRE(FILE='autoflat.fts',EXIST=EXIST)
       IF( EXIST ) THEN
         OPEN(26,FILE='autoflat.fts')
         CLOSE(26,STATUS='DELETE')
       ENDIF
C Vytvoreni vystupniho obrazku
       CALL FTINIT(26,'autoflat.fts',1,STATUS)
C zapsani informace do hlavicky
! dimenze a ostatni parametry se berou podle posledniho obrazku.
! pokud maji obrazky pretim odlisne dimenze, vysledek je nedefinovan 
       EXTEND = .FALSE.
       CALL FTPHPS(26,BITPIX,NAXIS,NAXES,STATUS)
       CALL FTPKYS(26,'OBJECT','FLAT-FIELD',' ',STATUS)
       CALL FTPKYE(26,'BSCALE',BSCALE,10,' ',STATUS)
       CALL FTPKYE(26,'BZERO',BZERO,10,' ',STATUS)
       CALL FTPKYE(26,'SKYMOD',SKYMOD,10,'Mean level',STATUS)
! komentar
       CALL FTPCOM(26,' ',STATUS)
       CALL FTPCOM(26,'Carefully!',STATUS)
       BUF='This file has been computed of the experimental code:'
       CALL FTPCOM(26,BUF,STATUS)
       CALL FTPCOM(26,VERSION,STATUS)

       print *,bitpix,naxis,naxes,status
C zapsani  dat
       CALL FTP2DJ(26,1,MAXX,NAXES(1),NAXES(2),IMAGE,STATUS)
C Zavreni souboru
       CALL FTCLOS(26,STATUS)
       IF( STATUS.NE. 0 )THEN
         PRINT *,'STATUS103=',STATUS
         STOP
       ENDIF

       print *,'Bye.'
       END       



